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Numerical Simulation of Three-Dimensional Self- Gravitating Flow 


John V. Shebalin* 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23681, USA 

The three-dimensional flow of a self-gravitating fluid is numerically simulated using a Fourier pseu- 
dospectral method with a logarithmic variable formulation. Two cases with zero total angular momentum 
are studied in detail, a 32 3 simulation (Run A) and a 64 3 simulation (Run B). Other than the grid size, the 
primary differences between the two cases are that Run A modelled atomic hydrogen and had considerably 
more compressible motion initially than Run B, which modelled molecular hydrogen. The numerical results 
indicate that gravitational collapse can proceed in a variety of ways. In the Run A, collapse led to an 
elongated tube-like structure, while in the Run B, collapse led to a flatter, disk-like structure. 


* Research supported by the National Aeronautics and Space Administration. This work was done while the author was in 
residence as a Visiting Scientist at the Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley 
Research Center, Hampton, VA 23681. 
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1 Introduction 


In recent years, Fourier methods have been used to simulate collapse of a gas due to a newtonian gravitational 
self-interaction, both in one-dimension (ID) and in two-dimensions (2D). In the ID case [1], the flow was 
assumed to be isentropic, 64 grid points were used, and it was found necessary to modify the values of density 
at each time-step by adding a spatial constant so as to ensure sufficient positive-definiteness of the density 
at all grid points. In the 2D case [2], the flow was again assumed isentropic, the maximum grid size was 
256 2 , the density was again ‘modified’ to ensure numerical stability, and, in addition, the flow was initially 
turbulent (with the notable result of hinderng collapse). 

Here, this work is extended in four ways. First, the full polytropic gas-dynamic equations are solved 
with no isentropic approximation. (Collapse is assumed to take place on length and time scales such that 
no overall cosmological expansion need be imposed on the dynamics.) Second, simulations are performed on 
three-dimensional grids (32 3 and 64 3 ). Third, the governing equations are written so that the logarithms of 
the density and temperature are solved for, rather than the density and temperature directly. This allows 
explicitly for the preservation of positive-definiteness of density, temperature, and pressure, since these are 
arrived at by exponentiation. Fourth, bulk viscosity is utilized to model more accurately the properties of 
molecular hydrogen. 

These extensions allow for more realistic simulations of the self-gravitational collapse of clouds of either 
atomic (H) or molecular {H 2 ) hydrogen, at least up to the point where density gradients are so steep that 
the computer code loses accuracy. In particular, two cases with zero total angular momentum are studied 
in detail, one on a 32 3 grid (Run A: H) and the other on a 64 3 grid (Run B: H 2 ). The primary difference 
between the two cases (other than grid size and chemical constituents) lies in the ratio of compressible kinetic 
energy to total kinetic energy. Run A had an initial ratio of 0.5, while Run B had an initial ratio of 0.1. 
Collapse from a constant density initial state is clearly different for the two cases: in Run A, the collapse is 
to a tube-like structure, while in Run B, it is to a flat, disk-like structure (a Zeldovich ‘pancake’ [3]). 

Following brief descriptions of the governing equations and numerical method, numerical results will 
be discussed. Flow parameters will be defined, as well as measures of anisotropy which are useful for 
quantitatively describing self-gravitational collapse. Finally, a summary and conclusion will be presented. 
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2 Physics &; Numerics 


Placing p - e x and T = t c into the standard polytropic gas equations yields the basic non-dimensional 
equations in a logarithmic formulation [4]: 
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Here, tj, = 5,u ; + <9 ; u, - 2/3 fyV • u and I = [fy] is the unit dyadic. In (2), $ is the newtonian gravitational 
potential which satisfies V 2 $ = k)p, where kj is the Jeans wavenumber [5]; for Run A, kj = y/2 and for 
Run B, kj = 1.5. The adiabatic index 7 is constant, and since the gas is either H or H 2 , 7 = 5/3 for Run A 
and 7 = 7/5 for Run B. The dimensionless shear viscosity ft and thermal conductivity k are also constant; 
the Prandtl number, Pr = 7/1/fC, has a value of ~ 1, which will determine k: k = 7^ (p = 0.02 for Run A 
and p. = 0.005 for Run B). The ratio of bulk to shear viscosity, /? = pb/p, will also be constant: 0 = 1 for 
Run A (to represent a small admixture of H 2 with primarily H) and 3 = 32 for Run B (only H 2 ) [6]. 

Details of the numerical method used here have been given previously [4, 7]. The only modification 
lies in the addition of the self-gravitation term -V$ in (2). Using V 2 $ = k 2 jp and Fourier transformation 
gives -V4> — ► ik 2 kk~ 2 p(k), which is very straightforward to implement when (1)— (3) are written in terms 
of Fourier coefficients (as is done in the method utilized here). We now proceed on to the results of the 
simulations. 


3 Results 

Run A used a 32 3 spatial grid and ran a total of 15300 time-steps (A f's), with 2 c.pu-sec/A t on a Cray AMP. 
Run B used a 64 3 spatial grid and ran a total of 49800 A t's, with 13 cpu-sec/A t y also on a Cray YMP. 
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Run A was terminated at simulation time t = 6.57 when the dissipation wave number kr> (as defined in [4]) 
had increased so as to equal the maximum wave number k max = 16. Similarly, Run B was terminated at 
simulation time t = 5.01 when the dissipation wave number ko had increased so as to equal its maximum 
wave number k ma x = 32. The microscale Reynolds number R\ (also defined in [4]) ranged between 17 and 
32 for Run A, and between 90 and 170 for Run B. 

In the following, ( Q ) will denote the volume average (or average value per grid point) of any quantity Q. 
Thus, the various energies associated with the collapsing clouds of gas are the kinetic energy Ek = (pu 2 / 2), 
the internal energy Ei - {pT/[y(y- 1)]), the gravitational energy Eg = {p$/2), and the total energy 
E — Ek + Ej + Eq> The time evolution of these energies is shown in Figure 1 for Runs A and B. 

As is evident in Figure 1, collapse did not begin immediately, but rather after the initial turbulent motion 
had decayed to a lower level, in accord with previous results [2]. (The initial turbulent velocity spectra was 
|ti(k )| 2 ~ k A exp — 2fc 2 /fc 2 , with k 0 = 3 for Run A and k 0 = 4 for Run A.) However, once collapse had 
progressed sufficiently, both the mean-squared dilatation $ = y^(V *u) 2 ^ and the mean-squared vorticity 
(the enstrophy) fi = |^(V x u) 2 ^) began to increase. This is shown in Figure 2, where it is clear that 
Run A had more compressible velocity initially than Run B did. [The compressible part of the velocity is 
u c (k) = fc“ 2 kk • u(k).] Note that both and Q begin to increase near the end of the runs. 

The relative amount of compressible velocity in a flow is quantified by the ratio \ = (u?) / (u 2 ), whose 
importance was first recognized by Passot and Pouquet [8]. The time evolution of \ for Runs A and B 
is shown in Fiqure 3. For Run A, x — 0*25 before collapse initiates, while for Run B, x — 0-02 during 
the equivalent time period. Although there are many differences between Runs A and B, the qualitatively 
different behavior of x prior to collapse may be related to the qualitatively different dynamical behavior of 
the gas in Runs A and B during the collapse itself. 

The results of gravitational collapse in Runs A and B are seen in Figures 4 and 5. These figures represent 
the densities of the respective clouds of gas at the end of each run and are drawn by finding the point of 
maximum density in either run, extracting values of density on the three orthogonal planes (x — y, y — z, and 
z - x) which contain this point, and then shifting the points on these periodic planes so that the maximum 
density lies in their centers. In Figure 4 it is clear that the gas has collapsed to a tube-like structure; in 
Figure 5 it is clear that the gas has collapsed to a much flatter, somewhat irregular, disk-like structure. 
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In order to get a better picture of how these structures developed, it is useful to define the following 
measures of anisotropy Mj (j — x,y,z ) [9]: 


Mj = 


3 (m 3 ) i 

2((Vp) 2 ) 2 


( 4 ) 


The time evolution, of the Mj for Runs A and B are shown in Figure 6. Note that the differences in the final 
states of Runs A and B, as shown in Figures 4 and 5, grew directly from their initial states, rather than 
going through a qualitatively different intermediate stage. 


4 Summary & Conclusion 

In this work a three-dimensional pseudospectral computer code was used to simulate the initial stages of 
the gravitational collapse of clouds of either atomic or molecular hydrogen. The full polytropic gasdynamic 
equations were solved, with the primary dynamical variables being In/?, lnT, and u. In addition to shear 
viscosity, bulk viscosity was also utilized to more exactly model the properties of molecular hydrogen. 

Two runs were examined in detail, Run A (on a 32 3 grid) and Run B (on a 64 3 grid). As the figures 
show, the nature of the collapse was different for each run: in Run A, the collapse proceeded to a quasi- ID 
structure, while in Run B, the collapse proceeded to a quasi-2D structure. One primary difference between 
the runs lay in the relative amount of compressible turbulent kinetic energy prior to collapse. This suggests 
that turbulent density fluctuations may have an important role to play during gravitational collapse. 

There is a large amount left unexplored in all of this, with respect to variation in initial conditions, 
transport properties, and grid size, or with respect to the effects of magnetic fields or radiation, for example. 
The sheer length of a single, moderately resolved simulation (64 3 — » ~ 200 Cray YMP cpu-hours) precludes 
any quick but thorough investigation of ‘parameter space.’ Instead, it is clear that only a slow and patient 
effort (along with exponential increases in computer power) will yield the results and insight which are 
sought. 
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Figure 2. Enstrophy Q and mean-squared dilatation for Runs A and B. 



Figure 3. Compressibility ratio % for Runs A and B. 
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Figure 5. Density contours in orthogonal planes intersecting 

the point of maximum density ( = 0.79) at t=5.01 

for Run B. 
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